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Abstract. 


Symplectic  Runge-Kutta  schemes  for  integration  of  general  Hamiltonian  systems  are 
implicit.  In  practice  the  implicit  equations  are  often  approximately  solved  based  on 
the  Contraction  Mapping  Principle,  in  which  case  the  resulting  integration  scheme 
is  no  longer  symplectic.  In  this  note  we  prove  that,  under  suitable  conditions,  the 
integration  scheme  based  on  an  n-step  successive  approximation  is  0(Sn+2)  away  from 
a  symplectic  scheme  with  <5  €  (0, 1).  Therefore,  this  scheme  is  “almost”  symplectic 
when  n  is  large. 

AMS  subject  classification:  65L05,  65L06,  65P10. 

Key  words:  Hamiltonian  systems,  symplectic  Runge-Kutta  methods,  Contraction 
Mapping  Principle,  successive  approximation. 


1  Introduction. 


Geometric  integration  methods,  numerical  methods  that  preserve  geometric 
properties  of  the  flow  of  a  differential  equation,  outperform  the  off-the-shelf 
schemes  (e.g.,  fourth  order  explicit  Runge-Kutta  method)  in  predicting  the  long¬ 
term  qualitative  behaviors  of  the  original  system  [5] .  An  important  class  of  geo¬ 
metric  integrators  are  symplectic  integration  methods  for  Hamiltonian  systems 
[9] .  Consider  a  Hamiltonian  system 


(1.1) 


m  =  - ^ 
m  = dJ^  ’ 


with  the  Hamiltonian  H(p,q),  where  (p,q)  G  x  Q  for  some  integer  d  >  1, 
and  Q,  the  configuration  space,  is  some  d-dimensional  manifold.  For  ease  of 
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discussion,  in  this  note  we  assume  Q  =  Rd,  but  the  results  we  present  here 
apply  to  the  case  of  a  general  Q  directly.  Let  2  =  (p,  q ),  the  system  (1.1)  can  be 
rewritten  as: 


(1.2) 

where 


A 


z(t )  =  =  JS7zH{z{t)) 


J  = 


0  -Id 
Id  0 


Id  denotes  the  d-dimensional  identity  matrix,  and  V2  stands  for  the  gradient 
with  respect  to  z. 

When  the  Hamiltonian  has  a  seperable  structure,  i.e.,  H(q,p )  =  T(p )  +  V{q), 
explicit  Runge-Kutta  type  algorithms  exist  which  preserve  the  symplectic  struc¬ 
ture  [4,  11,  3,  7].  However,  this  is  not  the  case  for  general  Hamiltonian  systems. 

An  s-stage  Runge-Kutta  method  to  integrate  (1.2)  is  as  follows  [6]: 


(1.3) 


Vi  =  z0  +  t  aijfiVj),  i=  1,  • 
zi  =  zo  +  r^jWfe) 


',s 


where  r  is  the  time  step,  Zo  is  the  initial  value  at  time  1 0,  z\  is  the  numerical 
solution  at  time  to  +  r,  an] ,  bi  are  appropriate  coefficients  satisfying  the  order 
conditions  of  the  Runge-Kutta  method. 

Let  vIfT  be  the  one  time-step  flow  associated  with  the  algorithm  (1.3),  i.e., 
Zi  =  ’Lr(zo)-  From  [8],  the  transformation  'hr  preserves  the  symplecticness  of 
the  original  system  (1.2)  if 


(1.4) 


bi&ij  -I-  bjCtji  bibj  —  0,  i ,  j  —  1, 


Thus  if  (1.4)  is  satisfied,  we  have: 
(1.5) 


=  0, 

oz0  dz0 


where  T  denotes  the  transpose.  The  condition  (1.4)  forces  the  symplectic  Runge- 
Kutta  method  (1.3)  to  be  implicit.  In  the  interest  of  computation  efficiency, 
Aubry  and  Chartier  investigated  pseudo-symplectic  Runge-Kutta  methods,  which 
are  explicit  and  conserve  the  symplectic  structure  to  a  certain  order  [1] .  We  also 
note  the  closely  related  work  in  [2] ,  where  the  error  estimate  for  the  Lie-Poisson 
structure  is  established  for  integration  of  Lie-Poisson  systems  using  the  mid¬ 
point  rule. 
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In  this  note,  we  take  a  different  approach  from  [1],  Successive  approximation 
based  upon  the  Contraction  Mapping  Principle  is  often  used  to  obtain  an  ap¬ 
proximate  solution  to  yi  in  (1.3).  The  resulting  integration  scheme  based  on  the 
approximation  is  no  longer  symplectic.  It’s  of  interest  to  investigate,  to  what 
extent,  the  symplectic  structure  (1.5)  is  preserved  by  the  approximation  scheme. 
The  rest  of  this  note  is  devoted  to  answering  this  question,  and  it  turns  out  that 
the  scheme  using  an  n-step  approximation  is  0(Sn+2)  away  from  a  symplectic 
one  with  0  <  S  <  1.  Therefore,  when  n  is  large  enough,  the  approximation 
scheme  is  “almost”  symplectic. 


2  A  successive  approximation  method. 


Denote 


A 

y  = 


1  yi 

\  Vs 


F(y)  = 


1  f(yi) 

\  f(Vs) 


b  =  (6i ,  •  •  • ,  ba),  A0  =  [ay-],  and  A  =  Ao^hd,  where  “(g)”  denotes  the  Kronecker 
(tensor)  product.  We  recall  for  two  matrices  M  =  [to,, ]  and  R  =  [r,  ,■ ] ,  the 
Kronecker  product 


M  g>  R  = 


muR  777.12-R 
TO21-R  77122-R 


The  algorithm  (1.3)  can  now  be  written  as 


y  =  G(z0,  y)  =  1  <g>  z0  +  rAF(y) 
Z!  =  20  +  rb  <g>  /2dF(y) 


where  1  is  an  s-dimensional  column  vector  with  1  in  every  entry. 

As  noted  in  Section  1,  when  (1.4)  is  satisfied,  the  first  equation  in  (2.1)  is 
implicit  for  each  fixed  Zo ■  One  algorithm  often  used  to  solve  implicit  equations, 
is  the  successive  approximation  scheme  based  on  the  Contraction  Mapping  Prin¬ 
ciple  (see,  e.g.,  [10]): 

Lemma  2.1  (Contraction  Mapping  Principle).  Let  S  be  a  closed  subset 
of  a  Banach  space  X  and  let  p  be  a  mapping  that  maps  S  into  S.  If3p  £  (0, 1), 
such  that 


ll<p(£)  -  <p(y) II  <  p\\x-y\\,\/x,y  £  S, 
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then 


1.  there  exists  a  unique  x*  £  S  satisfying  x*  =  <p(x*); 

2.  x*  can  he  obtained,  by  the  method  of  successive  approximation 

a>+1l  =  ^(z[n]), 

starting  from  an  arbitrary  in  S;  and 

3.  the  approximation  error  satisfies  \\x^  —  x*||  <  —  a:*||. 

In  this  note  we  will  use  ||  •  ||  to  denote  the  norm  (or  the  induced  norm)  of  a 
vector,  matrix,  or  high  order  tensors,  and  the  precise  meaning  should  be  clear 
from  the  context.  The  following  proposition  shows  that  when  the  step  size  r  is 
small  enough,  for  each  fixed  zo,  the  first  equation  in  (2.1)  has  a  unique  solution 

y*: 

Proposition  2.2.  Let  O  C  R2d  be  a  bounded  open  set.  Let  f  be  locally 
Lipschitz  continuous.  Then  for  any  S  £  (0,1),  e  >  0,  there  exists  r(fl,e,S)  >  0 
dependent  on  0,e  and  S,  such  that,  Vr  <  r(f2,  e,  <5),  Vzo  £ 

1.  there  exists  a  unique  solution  y*  =  y*(zo)  for  the  first  equation  in  (2.1); 


2.  y*  can  be  approximated  by  successive  approximation 

f  y[nI  =  G(z0,y[n_1l)  . 

^  |  yl°]  =  1  <g>  Z0 

and 


3.  ||y["]-y*||<J"||y[°]-y*||. 

Proof.  Denote  N(Q,e)  the  e— neighbourhood  of  fi,  defined  as 
N(Cl,  e)  =  {z  £  K2d  :  min  \\z  —  zo||  <  e}, 

Zq  G 

where  denotes  the  closure  of  fl.  Denote  Ns(Q,e)  the  product  of  s  copies  of 
N(Ll,e),  i.e., 

Ns(n,  e)  =  N(Sl,  e)  x  •  •  •  x  N(fl,  e). 

Since  N(Q,  e)  is  compact,  /  is  bounded  and  Lipschitz  continuous  with  some 
Lipschitz  constant  Lf  on  N(Cl,e).  Thus  there  exists  ri  >  0,  such  that  when 
t  <  n,  for  each  fixed  z o  £  O,  G(zq,  •)  maps  Ns(Ll,  e)  into  itself. 


5 


For  any  Zo  €  f 2,  for  y,  y'  €  Ns(tt,e),  by  the  definition  of  G, 

l  f(yi)  -  f(y'i) 

||G(^0,y)  -  G(z0,y')ll  =  IM  i 

\  f{ys)-f{y's) 

<  rLf\\A II  lly  —  y'll- 


For  <5  €  (0,1),  let  r2  =  xjpjj-  Now  for  r  <  r(fl,e,S)  =  min{ri,r2},  G(_20,-)  is 
a  contraction  mapping  for  each  fixed  zq  £  O.  All  the  claims  then  follow  from 
Lemma  2.1.  Note  that  r(Q,  e,  (5)  depends  on  Q,  e  and  <5.  □ 

Similarly  we  can  prove: 

Proposition  2.3.  Let  /  &e  globally  bounded  and  Lipschitz  continuous.  Then 
for  any  6  €  (0, 1),  there  exists  t(S)  >  0  dependent  on  5  only,  such  that,  Vr  < 
t(5),  for  each  fixed  zo  £  R2d,  G(2o,')  is  a  contraction  mapping  and  the  claims 
in  Proposition  2.2  hold. 

Remark  2.1.  Traditionally  implicit  Runge-Kutta  methods  have  been  used 
mostly  for  stiff  problems,  where  the  Lipschitz  constant  for  f  is  relatively  large  and 
the  convergence  of  successive  approximation  based  on  the  Contraction  Mapping 
Principle  is  slow.  However,  in  the  new  context  of  symplectic  integration,  we  are 
dealing  with  implicit  methods  even  for  nonstiff  problems.  Hence  the  successive 
approximation  plays  an  important  role  in  solving  the  implicit  equations. 

As  we  see  from  Proposition  2.2,  when  r  is  sufficiently  small,  the  solution  y* 
to  the  first  equation  in  (2.1)  is  a  function  of  z o,  and  we  can  write  it  as  y*(zo)- 
If  /  is  differentiable,  we  have  from  the  Implicit  Function  Theorem  that 

(2.3)  -jM^b)  =  [hsd  -  tA—  (y*^))]-1^  ®  I2d). 

azo  o  y 


3  Main  result. 


An  explicit  but  approximate  algorithm  to  solve  (2.1)  is  as  follows:  for  some 
n  >  0, 

f  y[fc]  =  G(2o,y[fc-1]),  k=  1,  •  •  • ,  n 

(3.1)  <  y[°l  =  1  (g)  20 

1  4"1  =  ^o  +  rb®  J2dF(yN) 

Remark  3.1.  The  scheme  (3.1)  based  on  n—step  successive  approximation  (to 
y* )  is  essentially  an  s(n  +  l)-stage  explicit  Runge-Kutta  scheme  with  coefficients 
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A  and  b,  where 


0 

1  0 


®  A>,  b=  (0,  •••,0 M,"*->ba). 


1  0 


Note  that  in  (2.1)  and  (3.1),  y*,  {y^}fc_0>  z\  and  (and  smooth  functions 
of  them)  are  all  continuously  differentiable  functions  of  Zq  if  /  is  differentiable 
and  t  is  sufficiently  small.  In  the  sequel  when  we  write,  e.g.,  or  ^■F(y[nl), 
we  think  of  y*  or  F(y["l)  as  a  function  of  zo  although  it  is  not  explicitly  written 
out. 

Denote  by  \E^  the  one  time-step  flow  associated  with  the  algorithm  (3.1),  i.e., 
(zq).  We  now  want  to  study  how  far  is  away  from  a  symplectic 
transformation.  The  following  lemma  will  be  essential  in  the  proof  of  our  main 
result  Theorem  3.3: 

Lemma  3.1.  Let  f l  c  K2d  be  bounded,  convex  and  open.  For  e  >  0,  let 
N(Q,e)  be  the  e-neighbourhood  of  Ll,  as  defined  in  the  proof  of  Proposition  2.2. 
Assume  that  f  is  twice  continuously  differentiable  on  Then  for  any 

6  £  (0,1),  there  exists  r(Cl,e,6)  >  0  dependent  on  fi,e  and  S,  such  that  when 
t  <  r(fl,e,(5),  for  each  fixed  z o  €  fi,  the  first  equations  in  (2.1)  and  (3.1)  have 
(unique)  solutions  y*  £  Ns(fl,e )  and  y^"l  £  Ns(fl,e),  respectively;  and 


(3.2) 


<  C^.e)**"*1, 


(3.3)  ||^-(F(yW)-F(y*))||<C,(f2,e)^+1,, 

oz  0 

where  C(f2,  e),  C'ffl,  e)  >  0  are  constants  dependent  only  on  f l  and  e. 

PROOF.  Since  /  is  differentiable,  it  is  Lipschitz  continuous  on  the  convex  set 
IV(fi,e).  By  Proposition  2.2,  there  exists  ri(fl,e,(5)  >  0,  such  that  when  r  < 
ri(fl,  e,  (5),  for  any  zo  £  fi,  G(zo,  •)  is  a  contraction  mapping,  y*,  y^  £  Ns(Ll,  e), 
Vfc  >  0,  and  (recall  (2.3)) 

(3-4)  ll^ll<C'1(Q,e), 

OZo 


^y[0]  <9y*  <9F  <9y* 

1^-  -  ^-||  =  \\rA^-(y*)^-\\  <  rC2( fl,e), 
dz0  dz0  dy  dz0 


(3.5) 
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where  C,( Cl,  e)  >  0,  *  =  1,  2,  are  constants  dependent  on  f l,  e. 
From  (2.1)  and  (3.1), 


(3.6) 


yN  _  y*  =  T^4(F(y["_1l)  —  F(y*)). 


Taking  derivative  of  both  sides  of  (3.6)  with  respect  to  Zo  and  re-arranging  terms, 
we  get 

dy ["1  dy* 
dz0  dz0 


(3.7) 


=  tA[ 


dy 


dzn 


dy 


dy 


dzn 


Denoting 


0M  £  rA^-(yW),  r[fc]  =  r7l(^(y[fc])  -  f-(y*))^, 

ay  ay  ay  ozq 


we  derive  from  (3.7) 


n—  1  n—1 


u  u  k=0  u  u  k—0  i=k+ 1 

which  implies 

k—0  k—0  i=k+l 

The  following  two  observations  are  in  order: 

1. 

(3.9)  ||0[fcl||<r||A||  max  ||  — (y)||  <  rC3(n,  e), 

yeAT»(n,e)  dy 

where  Cs(Cl,  e)  >  0  is  a  constant  dependent  only  on  Q  and  e. 
2.  When  r  <  ri(f2,  e,  (5), 

4T7' 

r\k] 


(3.10)  ||_(yW)-^(y*)||  <  max  J  ^-(y)||  ||y- -  y 


dF , 

dy ' 


d2F , 

yeiv»Tn,e) 1  dy2 
<  C'4(f2,e)5fc||y[°]-y*||, 


where  C^(il.  e)  >  0  is  a  constant  dependent  only  on  f l  and  e.  Combining 
(3.4)  and  (3.10),  and  using 

||y[0]  —  y* ||  <  r||A||  max  ||F(y)||, 

yeAfs(n,e) 


(3.11) 


we  have 

(3.12)  nr!*]  ||  <  T26kC5(Cl,e), 

for  some  constant  Cffil,  e)  >  0. 


Pluggin  (3.5),  (3.9)  and  (3.12)  into  (3.8),  we  obtain  after  some  manipulations 

^^iTC^,e){rC3Mr+ 


A 


We  now  let  T2(fi,  e,  S)  =  2c3(0  7)  ’  anc^  ^ 

r( ft,  e,  (5)  =  min{n(f2,  e,  <5),  r2(fi,  e,  (5)}. 


It’s  easy  to  verify  that,  Vr  <  r(0,  e,  <5), 


II 


9yN 

9^0 


<  C,(n,e)(Jn+1, 


where  C(Si,e)  =  $$1%  +  2cf(n%  •  This  Proves  (3-2)- 
To  show  (3.3),  we  note  that 


/-(F(yH)  -  F(y*))  =  5(yI»l)(|H  -  §1)  +  (§V)  -  §V))§^. 
ozq  dy  ozo  ozo  dy  ay  dzo 

and  then  use  (3.2),  (3.4),  (3.10)  and  (3.11).  □ 

Similarly  we  can  prove: 

Lemma  3.2.  Let  f  be  globally  bounded  and  twice  continuously  differentiable, 
with  bounded  first  order  and  second  order  derivatives.  Then  for  any  S  €  (0, 1), 
there  exists  t(6)  >  0  dependent  on  6  only,  such  that  when  r  <  t(S),  for  each 
fixed  zq  £  R2  ,  the  first  equations  in  (2.1)  and  (3.1)  have  (unique)  solutions  y* 
and  yM,  respectively;  and 


(3.13) 


9ylnl  dy* 


dzn 


dzn 


<  CSn+1, 


(3.14)  ||^-(F(yN)  -F(y*))||  <  C'Sn+1, 

dzo 

for  some  constants  C,  C'  >  0. 

We  are  now  ready  to  present  the  main  result  of  this  note: 

Theorem  3.3.  Let  O  c  R2d  be  bounded,  convex  and  open.  For  e  >  0, 
let  N(fl,e)  be  the  e-neighbourhood  of  Ll.  Assume  that  f  is  twice  continuously 
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differentiable  on  N(fl,e).  Consider  the  algorithm  (3.1)  and  let  \E'i^  be  the  one 
time-step  flow  associated  with  (3.1).  Let  (1-4)  be  satisfied.  Then  for  any  S  € 
(0,1),  there  exists  r(fl,e,S)  >  0  dependent  on  Q,e  and  5,  such  that  when  r  < 
r(Ll,e,S), 


(3.15)  ||( 1  (zo)  )T  J( ^ 1  (~o) )  -  J||  <  C(Q,  e)5n+2 ,  Vz0  £  O, 

Oz  0  OZo 

where  <7(11,  e)  is  a  constant  dependent  on  O  and  e. 

PROOF.  By  Lemma  3.1,  we  can  find  r(Cl,e,S)  >  0,  such  that  when  r  < 
r(Q,e,i 5), 

(3.16) 


d 


||  (F(y  W)  _  F(y*))||  <  e)Sn+1, 

OZo 

for  some  constant  Ci(f2,  e),  where  y*  and  yW  are  solutions  to  the  first  equations 
in  (2.1)  and  (3.1),  respectively.  Let  'hr  be  the  one  time-step  flow  associated  with 
(2.1).  From  (2.1)  and  (3.1),  we  have 


AH(~o)  =  M(*o)  -  *T(z0)  =  rb  ®  /2d(F(yN)  -  F(y*)), 


which,  by  (3.16)  and  the  definition  of  r(fi,  e,  d) 


(3.17) 


SAW^o), 


dzn 


<C2(ft,e)Sn+2,  Vz0  e  fi, 


where  the  constant  C2(fb  e)  depends  only  on  Q  and  e.  We  now  write 

dA^o)  |  d^T(z0)^T  dAM(z0)  |  G>frr(zo)) 

dzn  dzt)  dzn  dzn 


< 


,dAW(zo)  T  dAWjzp) 
dzn 


,dAW  (go)  r  T,d*r(zo). 
[  dzo  >  [  dzo  ' 


dz0  oz0 

d^r(z0)^T  9AN(£o)  ..gjbKfo)  r  gjM£o) 

dz0  )  [  dzo  jll  +  IK  dzo  ’  [  dzo  ’ 


where  the  last  term  vanishes  when  (1.4)  is  satisfied. 
Finally,  we  note  that 


(3.19)  ||^^)  ||  =  ||/2d  +  rb®/2d^(y*)^||  <  C3(fi,e) 

dzo  dy  dzo 

for  some  constant  Cs(Ll,e)  >  0,  where  (3.4)  is  used.  Combining  (3.17),  (3.18), 
and  (3.19)  yields  (3.15).  □ 
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A  global  version  of  Theorem  3.3  can  be  proved  analogously: 

Theorem  3.4.  Let  f  be  bounded  and  twice  continuously  differentiable,  with 
bounded  first  order  and  second  order  derivatives.  Consider  the  algorithm  (3.1) 
and  let  be  the  one  time-step  flow  associated  with  (3.1).  Let  (1-4)  be  satisfied. 
Then  for  any  5  G  (0,1),  there  exists  t(6)  >  0  dependent  only  on  S,  such  that 
when  t  <  r(<5), 


(3.20) 


dz0 

for  some  constant  C  >  0. 


^1(*o))TJ(d*f(2o))  -  J||  <  C6n+2,  Vzo  G  R2d, 


dz0 
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